package mylib.dir;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.Random;

/**
 * Created with IntelliJ IDEA.
 * User: 1
 * Date: 16.07.12
 * Time: 18:31
 * To change this template use File | Settings | File Templates.
 */
public class MathUtils {
    private static Random random = new Random(23095823597L);

    public static boolean[] eratosthenes(int up) {
        boolean[] isPrime = new boolean[up + 1];
        Arrays.fill(isPrime, true);
        isPrime[0] = isPrime[1] = false;
        for (int i = 2; i * i <= up; ++i) {
            if (isPrime[i]) {
                for (int j = i + i; j <= up; j += i) {
                    isPrime[j] = false;
                }
            }
        }
        return isPrime;
    }

    public static long trivialEulerFunction(long x) {
        long res = x;
        for (long i = 2; i * i <= x; ++i) {
            if (x % i == 0) {
                while (x % i == 0)
                    x /= i;
                res -= res / i;
            }
        }
        if (x != 1)
            res -= res / x;
        return res;
    }

    public static long fastEulerFunction(long x) {
        if (x <= 1)
            return 0;
        Pair<Long, Integer>[] f = factorize(x);
        long r = 1;
        for (Pair<Long, Integer> pair : f) {
            long cur = 1;
            for (int j = 0; j < pair.getSecond(); ++j)
                cur *= pair.getFirst();
            r *= (cur - cur / pair.getFirst());
        }
        return r;
    }

    public static int[] getPrimesWithSieve(int up) {
        boolean[] isPrime = eratosthenes(up);
        ArrayList<Integer> arrayList = new ArrayList<Integer>();
        for (int i = 2; i <= up; ++i)
            if (isPrime[i])
                arrayList.add(i);
        return ArrayUtils.toArray(arrayList);
    }

    public static int[][] binomial(int len, int modulo) {
        int[][] cnk = new int[len + 1][len + 1];
        for (int i = 0; i < cnk.length; ++i) {
            cnk[i][i] = cnk[i][0] = 1;
            for (int j = 1; j < i; ++j)
                cnk[i][j] = (cnk[i - 1][j] + cnk[i - 1][j - 1]) % modulo;
        }
        return cnk;
    }

    public static int[][] binomial(int len) {
        return binomial(len, Integer.MAX_VALUE);
    }

    public static int[] factorials(int len, int modulo) {
        int[] f = new int[len + 1];
        f[0] = 1 % modulo;
        for (int i = 1; i < f.length; ++i) {
            f[i] = (int) ((long) f[i - 1] * i % modulo);
        }
        return f;
    }

    public static int[] reverseFactorials(int len, int modulo) {
        int[] f = new int[len + 1];
        f[0] = 1;
        for (int i = 1; i < f.length; ++i) {
            if (MathUtils.gcd(i, modulo) != 1)
                throw new IllegalArgumentException("this is not field modulo = " + modulo);
            f[i] = (int) ((long) f[i - 1] * inverse(i, modulo) % modulo);
        }
        return f;
    }

    private static long[] primesStack = new long[70];
    private static int pointer = 0;

    public static Pair<Long, Integer>[] factorize(long x) {
        if (x <= 0)
            throw new IllegalArgumentException("!!! number must be positive!!!");
        pointer = 0;
        if (x == 1) {
            primesStack[pointer++] = x;
        } else {
            recursiveFactorization(x);
        }
        Arrays.sort(primesStack, 0, pointer);
        int different = 0;
        for (int i = 0, j; i < pointer; i = j) {
            for (j = i; j < pointer && primesStack[i] == primesStack[j]; ++j) ;
            ++different;
        }
        Pair<Long, Integer>[] res = new Pair[different];
        different = 0;
        for (int i = 0, j; i < pointer; i = j) {
            for (j = i; j < pointer && primesStack[i] == primesStack[j]; ++j) ;
            res[different++] = new Pair<Long, Integer>(primesStack[i], j - i);
        }
        return res;
    }

    private static void recursiveFactorization(long x) {
        while (x % 2 == 0) {
            primesStack[pointer++] = 2;
            x /= 2;
        }
        if (x == 1)
            return;
        if (fastPrime(x))
            primesStack[pointer++] = x;
        else {
            long divisor = rho(x);
            recursiveFactorization(divisor);
            recursiveFactorization(x / divisor);
        }

    }

    private static long rho(long x) {
        long g = x;
        while (g == x) {
            long u = rand63bits() % (x - 2) + 2;
            long v = u;
            long s = rand63bits() % (x - 2) % 2;
            g = 1;
            while (g == 1) {
                u = poly(u, s, x);
                v = poly(poly(v, s, x), s, x);
                g = gcd(u - v, x);
            }
        }
        return g;
    }

    private static long poly(long u, long s, long x) {
        return (handyMultiply(u, u, x) + s) % x;
    }

    public static boolean trivialPrime(long x) {
        if (x <= 1)
            return false;
        for (long i = 2; i * i <= x; ++i) {
            if (x % i == 0)
                return false;
        }
        return true;
    }

    public static boolean fastPrime(long n) {
        if (n == 2) return true;
        if (n <= 1 || n % 2 == 0) return false;
        long d = n - 1;
        int s = 0;
        while (d % 2 == 0) {
            d /= 2;
            s++;
        }
        for (int iter = 0; iter < 15; ++iter)
            if (!strictlyPseudoPrime(n, s, d))
                return false;
        return true;
    }

    private static boolean strictlyPseudoPrime(long n, int s, long d) {
        long a = rand63bits() % (n - 2) + 2;
        if (gcd(a, n) != 1) return false;
        a = binaryPowerHandyMultiply(a, d, n);
        for (int j = 0; j < s; ++j) {
            if (j == 0 && a == 1 || a == n - 1) return true;
            a = handyMultiply(a, a, n);
        }
        return false;
    }

    public static long binaryPowerHandyMultiply(long base, long power, long mod) {
        long res = 1 % mod;
        while (power > 0) {
            if (power % 2 == 1)
                res = handyMultiply(res, base, mod);
            base = handyMultiply(base, base, mod);
            power /= 2;
        }
        return res;
    }

    public static long handyMultiply(long a, long b, long mod) {
        a %= mod;
        b %= mod;
        long res = 0;
        while (b > 0) {
            if (b % 2 == 1) {
                res = res + a;
                if (res >= mod)
                    res -= mod;
            }
            a += a;
            if (a >= mod)
                a -= mod;
            b /= 2;
        }
        return res;
    }

    private static long rand63bits() {
        return (((long) random.nextInt(1 << 30)) << 31L) + (random.nextInt(1 << 30) << 1) + random.nextInt(2);
    }

    public static int binaryPower(int base, int power, int modulo) {
        int res = 1 % modulo;
        while (power > 0) {
            if (power % 2 == 1)
                res = (int) ((long) res * base % modulo);
            power /= 2;
            base = (int) ((long) base * base % modulo);
        }
        return res;
    }

    public static int gcd(int a, int b) {
        if (a < 0) a = -a;
        if (b < 0) b = -b;
        while (a > 0 && b > 0) {
            if (a >= b) a %= b;
            else b %= a;
        }
        return a + b;
    }

    public static long gcd(long a, long b) {
        if (a < 0) a = -a;
        if (b < 0) b = -b;
        while (a > 0 && b > 0) {
            if (a >= b) a %= b;
            else b %= a;
        }
        return a + b;
    }

    public static long[] gcdExtended(long a, long b) {
        if (a == 0) {
            return new long[]{b, 0, 1};
        } else {
            long[] tmp = gcdExtended(b % a, a);
            return new long[]{tmp[0],
                    tmp[2] - b / a * tmp[1], tmp[1]};
        }
    }

    public static int[] gcdExtended(int a, int b) {
        if (a == 0) {
            return new int[]{b, 0, 1};
        } else {
            int[] tmp = gcdExtended(b % a, a);
            return new int[]{tmp[0],
                    tmp[2] - b / a * tmp[1], tmp[1]};
        }
    }

    public static int inverse(int a, int m) {
        long x = gcdExtended((long) a, (long) m)[1];
        x = (x % m + m) % m;
        return (int) x;
    }

    public static long lcm(long a, long b) {
        return a / gcd(a, b) * b;
    }
}
